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ABSTRACT 

Planetary bodies form by accretion of smaller bodies. It has been suggested that a very efficient way to grow protoplanets is by ac- 
creting particles of size <Kkm (e.g., chondrules, boulders, or fragments of larger bodies) as they can be kept dynamically cold. We 
investigate the effects of gas drag on the impact radii and the accretion rates of these particles. As simplifying assumptions we restrict 
our analysis to 2D settings, a gas drag law linear in velocity, and a laminar disk characterized by a smooth (global) pressure gradient 
that causes particles to drift in radially. These approximations, however, enable us to cover an arbitrary large parameter space. The 
framework of the circularly restricted three body problem is used to numerically integrate particle trajectories and to derive their 
impact parameters. Three accretion modes can be distinguished: hyperbolic encounters, where the 2-body gravitational focusing en- 
hances the impact parameter; three-body encounters, where gas drag enhances the capture probability; and settling encounters, where 
particles settle towards the protoplanet. An analysis of the observed behavior is presented; and we provide a recipe to analytically 
calculate the impact radius, which confirms the numerical findings. We apply our results to the sweepup of fragments by a protoplanet 
at a distance of 5 AU. Accretion of debris on small protoplanets (<50 km) is found to be slow, because the fragments are distributed 
over a rather thick layer. However, the newly found settling mechanism, which is characterized by much larger impact radii, becomes 
relevant for protoplanets of ~10^ km in size and provides a much faster channel for growth. 

Key words. Planets and satellites: formation - Protoplanetary disks - Minor planets, asteroids: general 
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1. Introduction 

■ We consider how gas drag affects the collision rates between 
, a big body - a planetesimal or protoplanet - and small parti- 

■ cles, e.g., dust, chondrules, or bou l ders. Although the co re ac- 
I cretion model ( iPollacketal .11 19961: iHubickvi et al.ll2005h in its 
. initial stages, i.e., before the formation of a ~10 Earth mass (M®) 

core, concerns the accumulation of solid bodies, the role of the 
[ gas cannot be overstated. In the early phases of planet forma- 

■ tion - the growth of dust to planetesimals - the gas damps the 
' velocities of small particles. Initially, the (relative) velocities be- 
. tween particles are tiny and this is the reason why dust grains 
' can coagulate due to interm olecular forces (Dominik & Tielen^ 
I Il997t iBlum & Wurmll2000b - an effect much harder to envision 
, in the diffuse interstellar medium or even in molecular clouds. 

■ In this stage mechanisms that induce a relative velocity among 
' the dust particles include Brownian motion, settling, radial drift, 
1 and turbulent motions. The latter three are all functions of the 

particle's stopping time, a measure of how well particles couple 
to the gas. With increasing size (or, more correctly, increasing 
mass-to-surface area) particles couple less well to the gas and 
relative velocities increase, culminating in the so called meter- 
size barrier, which, at our current level of understanding, can 
best be overcome by the combined efforts of turbulent concen- 
tration an d gravitational collapse dJohansen et al.l l2007l l2009t 
ICuzzi et a l. 2010). 

Gas drag also affects the coUisional behavior at a much 
later stage, when protoplanets accrete planetesimals of perhaps 
~I - 10^ km in size. The collisional cross section between these 
big bodies is increased by gravitational focusing, i.e., the body 
can accrete particles at a cross section larger than its geometri- 



cal cross section due to gravitational deflection (ISafronovll 19691: 
IWetheriU & StewartllI989tlGreenzweig & Lissaueiill990ilI992l) 
This effect, however, is very sensitive to the velocity Ua at which 
the bodies approach: if y,, is too large, the focusing vanishes. In 
planetesimal accretion theory it is expected that a protoplanet 
will excite the random motions (eccentricities and inclinations) 
of the bodies it is accreting from, leading to a self-regulated ac- 
cretion behavior, which slows down the growth dlda & Makiii^ 
Il993t lKokubo&Idalll998t lOrmel et al.lHoiOah . Gas drag can 
provide some relief since, by damping the random motions of the 
planetesimals, the gravitational focusing is kept large. Moreover, 
the capture probability of planetesimals is also significantly 
increased when (protolplanets are surrounded by atmospheres 
(llnaba & Ikomal 120031: iTanigawa & OhtsukillloTob - again, gas 
drag is the mechanism that facilitates their accretion. Still, it is 
unclear if these effects are sufficient to overcome the timescale 
problem, i.e., to grow protoplanets t o ~10 Mm within the time 
the gas disk dissipates (^lO*" yr); see lLevison et al.l (1201 Ot) for a 
recent review. 



Due to the dynamical heating of planetesimals, planetesimal- 
planetesimal collisions may become di sruptive, producing 
smaller p lanetesimals or even fragments (IWetheriU & Stewart! 
119931 : iLeinhardt et al. 2009). These fragments can be kept dy- 
namically cold, e.g., by mutual collisions or by gas drag. The ac- 
cretion then takes place at low Vg - t he shear-dominated reg ime - 
which is very favorable for growth dGoIdreich et al.ll2004h . The 
generation of large amounts of fragments therefore can signif- 
icantly boost accretion . In particular, the accretion rate in the 
two dimensional (interactions are confined to a plane), gas-free, 
three-body regime (including the gravity of the central star) is 
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derived by a number of studies to be 



dM\ 



UafRHV,,^ 



(1) 



(e.g., 
I2OOI 



Ida & Nakazawa' '1989"; 'Greenbera et al.' 'l99Jt iRafikovl 
^AfeidenschiIling 2 005; Ormel et al. 2010b; the numerical 
constant is adopted from llnaba et al.ll200lh . where S is the den- 
sity in solids, Rf, the Hill radius. 



Ri, 



Mp_ 
3M* 



1/3 



(2) 



Vh the Hill velocity, y/, — Ri,0., a the semi-major axis, Q 
the corresponding orbital frequency, Mp/Mi, the ratio between 
the mass of the protoplanet and the central star, and Op 
the ratio between the protoplanet radius and the Hill radius, 
Op - Rp/Rfi- Equation ([T]) is often used in statistical models 
for the accretion rate (e.g.. Ilnaba etal]|2001 : IChamberi 20061: 



Brunini & Benvenutd I2OO8I: IChambersI 12008'; Kobavashi et al 
2010:). It represents a fast accretion rate. Kenyo n & Bromley 
(|2009|) . applying such a fragmentation-driven accretion scenario, 
calculate that the core formation process can be completed 
within 10* yr. 

How would gas drag affect these conclusions; i.e., does the 
rather large accretion rate of Eq. ([T]i also materialize in the pres- 
ence of gas drag? Qualitatively, two directions can be envi- 
sioned. On the one hand, the dissipative nature of the drag will 
enhance the collision (impact) radius, like in the case of a dense 
atmosphere. Conversely, strong particle-gas coupling will sup- 
press the accretion efficiency since the gas after all is not ac- 
creted but flows past the object (until the point where it has be- 
come more massive than 10 M^ and gas runaway accretion kicks 
in). It is a priori unclear which aspect of the drag - the coupling 
or the dissipation - will turn out to be the more important. 

To address these questions we include gas drag as an addi- 
tional force to the restricted 3-body problem that has been pre- 
viously used in calculating accretion rates in gas-free systems 
(or in systems where g as can be neglected; Petit & Henon 198^ 
llda & Nakazawairi989l) . Using appropriate scaling behavior, we 
show, in Sect. 12] that the system of equations containing all the 
physics can be restated into two dimensionless parameters: the 
dimensionless headwind velocity (^i that the protoplanet experi- 
ences and the dimensionless stopping time (Stokes number, St) 
of the particle. Our setup is idealized in the sense that we assume 
a steady gas flow of constant density (i.e., no pressure fluctua- 
tions or atmospheres), a drag law linear in velocity (applicable 
to small particles), and only consider drift motions of particles. 

After having outlined our setup in Sect.|2l Sect. [3] considers 
the geometrical limit, in which the 2 body interaction is absent 
or can be ignored. In Sect. |4] we perform an extensive parameter 
study to obtain the impact parameters as function of the relevant 
dimensionless quantities. Section|5] presents an analytic model 
to obtain the impact radii and accretion rates from first princi- 
ples, which we compare to our measured values. Section|6]illus- 
trates the significance of our result by calculating the protoplanet 
growth timescale in which we apply a correction to account for 
the scaleheight of the particle layer. We discuss limitations of 
our results and summarize in Sect. [7] 

2. Sketch of problem and approach 

2.1. Definition of impact radius 

In this study we will calculate both numerically and analytically 
impact radii, b^. In 3D systems, the collision rate dM/dt is the 



product of the velocity at which the bodies approach each other, 
the approach velocity, (;„, the cross section for collisions, o-, and 
the volume density in solids p that are accreted, dM/dt - pcrua- 
In 2D configurations the vertical dimension is lacking and we 
define 



idM\ 
['di'liD ~ 



Iba-Va'^ - Pco\^, 



(3) 



where Pcoi = '^.ba -Vg is the specific collision rate (cf. 
iNakazawa et al.|[l989l) . In the drag-free regime we indicate the 
impact radius b^ by^f. Although we primarily focus on 2D- 
configurations. Sect. |6l considers a 3D extension in which we 
apply the derived b^ also for the vertical dimension. 

In the gas-free regime particles enter the Hill sphere from 
orbits both interior and exterior to that of the planet, see Fig.[I] 
Therefore, b^i is associated with the lengthscale over which 
particles impact for one of these branches. However, particles 
can only enter at spec ific intervals, 1.7/?/, < \b\ < 2.5Ri, (e.g., 
iGreenberg et al.l Il99lh ; particles on impact parameters \b\ < 
1 .IRi, move on horseshoe orbits that do not enter the Hill sphere. 
The approach velocity Vg for the 3-body regime is defined as the 
average shear velocity {3bQ./2) over the above interval, i.e.. 



1 r^^MQfo 

Va = r-j do = i.Zvii. 

2.5Rh - l.lRi, Ji.vR,, 2 

Using Va - 3.2ui, and equating Eq. (O with Eq. (O gives 
ferf = IJal'-Ri, 



(4) 



(5) 



as the (effective) impact radius for accretion in the 2D gas-free 
regime. Note that since Op <s: 1, b„f <s: R^, which signifies that 
not every particle that enters the Hill sphere will collide. For this 
reason we distinguish between i>app = 2.5, the impact parameter 
at which particles approach (which is related to v„), and ba-, the 
impact radius that enters in the expression for the collision rate 
f col. The fact that b^- + Zjapp is peculiar to the three-body regime, 
where the gravity of the central star becomes important. 

2.2. The circularly restricted three body problem modified by 
gas drag 

We briefly review the circularly restrict ed three b ody problem 
using the framework of Hill's equations dHillll 18781) and include 
a drag term. The restricted three body problem assumes that the 
mass of the third body (M3) can be neglected with respect to 
the masses of the other two bodies, M\,M2. Furthermore, it is 
assumed that the orbits of the bodies are confined to a single 
plane and that these are circular for the two massive bodies. In 
our case the first body is the central star (M*), the second the 
(proto)planet {Mp), and the third the (test) particle m. We then 
consider the motion of the test particle in a coordinate system 
centered on and rotating with the motion of the planet. The re- 
sulting equations of motions for m in such a frame rotating with 
angular frequency Oo read 



dt 



^ = F - 2fio X V + d^r. 



(6) 



where r — {x, y, z) are the coordinates in the comoving frame, 
2Qo X V is the Coriolis acceleration and Q^r the centrifugal ac- 
celeration. The force per unit mass, F, acting on the third body 
consist of the solar gravity, Fsun = 0.^a, the 2-body force with 
the protoplanet. Fab = GMr/r^, and the drag force with the gas 
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Vfw 



1/gas = 3xQ/2 + vij^ 



Rh 





.Mp 



m 



Fig. 1. Sketch of particle trajectories in the comoving frame. We 
consider the motion of the third (test) particle m in the comoving 
frame of the second body {Mp, the planet) while including the 
gravity of the central star. In the gas-free limit zero-eccentricity 
particles can enter the Hill sphere from both the first and the 
third quadrant (black curves) but only from specific impact pa- 
rameters indicated by the hatched regions. Particles arriving at 
closer impact parameters move on horseshoe orbits. The mag- 
nitude and direction of the gas velocity Vgas as seen from the 
comoving frame is indicated by the dashed arrows. Particle tra- 
jectories including gas drag (solid gray arrows) can be anything 
depending on the properties of the particle and the gas. 

Fdiag. Expanding the Fjun term around oq enables us to linearize 
Eq. ^ to obtain 



It 



( movy + 3Qgx ' 



GM 



{ x\ 

y 

z , 



+ F 



drag ■ 



(7) 



Next, we rewrite Eq. (|7]l in dimensionless form by normalizing 
lengths to Hill radii /?/, (see Eq. (|2]i) and times to Qq'. The unit 
of velocity is then the Hill velocity, u/, = Rh^o- It can be shown 
that in Hill units GM - ^n^R^^ = 3. Dropping the z-term as we 
will treat planar configurations only, Eq. (|7]i, in Hill units, reads 



d\' 

IF 



2v' + 3x' 

y 



3x'/r'' , 



drag' 



(8) 



where F^^.^^ is related to Fdrag as F^^^^ = Fdrag//?,^i^o- 



2.3. The gas drag force 

The drag force, Fdiag, can be expressed in terms of a stopping 
time ts, 

AV (V - Vgas) 



•^ drag — — " 



(9) 



where Vgas is the velocity of the gas in the comoving frame and 
Av the velocity difference between that of the particle and the 



gas, see Fig.[T] Due to pressure support, the gas rotates slower 
than Keplerian by a magnitude Hhw - rjVK, where vk is the 
Keplerian velocity at disk radius a and rj a dimensionless quan- 
tity that gives the fracti onal deviation from the Keplerian motion 
dNakagawa et al.llI986h : 



^ = 



dPIda 
2aO?-pg 



Vk 



(10) 



with Cg the sound speed. In the comoving frame the headwind is 
directed towards negative y. However, we should correct for the 
Keplerian shear; thus. 



Vgas = (-Whw - -^X)ty, 



(11) 



where e^ is the unit vector in the y direction. 

For the drag force we consider several regimes. The stop- 
ping time for solid spheres of internal density p, for particles of 
increasing size s reads dWeidensc hilling 1977a): 



ts ^{ 



( PsS 

PgCg 
4p.,r 

Pg\y - VgasI 



(Epstein drag) 
(Stokes drag) 
(Quadratic drag) 



(12) 



where pg the density of the gas, and ^^fp the mean free path 
of the gas. For small particles the Epstein regime holds. The 
Stokes regime supersedes the Epstein regime for particle sizes 
s > 9Anfp/4. In both the Epstein and the Stokes regime the gas 
drag is linear with velocity and the stopping time reflects a par- 
ticle property. These are the regimes for which our study is ap- 
plicable. In the quadratic regime the stopping time becomes a 
function of the particle velocity since here Fdrag °^ |Av|^. In fact, 
there is a transition regime between the Stokes and quadratic 
drag regimes where stopping times are proportional to |Av|'^'^, 
which we have, for reasons of simplicity, ignored here (follow- 
ing JRafikov 2004) ■ 

As an (approximate) upper limit for |Av| we can take the 
headwind velocity, Uhw The transition between the Stokes and 
the quadratic drag regimes then occurs at a size of 



:90 m 



2yhw 
C, 



(13) 



/ gg \/ Pg \ / t'hw r 
\ 105 cm s-i / \ 10-10 g cm-3 / \ 30 m s-' / 



where we used /"mfp - 2x lO^^/p^ (in cgs units:^Nakaga wa et al.l 
|1986|) . Since we consider a drag law that is linear in velocity, our 
results are only applicable for particle sizes less than in,ax- In the 
inner disk (where the gas density is large) the results should be 
applicable to the sweepup of chondrule-like particles and m-size 
boulders. In the outer disk, pg is much lower and the particles for 
which our results are applicable include (small) planetesimals. 
Expressed in dimensionless form the drag law reads 



, ^ Fdrag 
drag ,. o 



1 
St 



-u'y-(u]-3x'/2 



(14) 



where we used Eq. (|9]l for Fdmg and Eq. ( fTTT ) for Vgas, normalized 
velocities to y/,, and have introduced the Stokes number, St = 
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Fig. 2. Relation between dimensionless and physical quantities. 
The Stokes number (St, solid lines) and the dimensionless head- 
wind velocity ^„, (dashed lines) are plotted on the (/-axis as func- 
tion of the radius s of the particle and the radius Rp of the 
protoplanet. Note the different units of s and Rp on the x-axis. 
The black dot denotes imax (Eq. (fT3])). Lines are shown for : (i) 
Cg - 10^ cm s"' and pg - 10"^ g cm"-' at a disk radius of 1 AU 
(black lines) and (ii) c^ = 6 x lO'* cm s"' andp^ - 10"" g cm""* 
at a position of 10 AU (gray lines). The internal density of solids 
is fixed at ps = 3 g cm"-', the nebula headwind is iihw = 30 m s"', 
and the mass of the central object is solar. 



Table 1. Dimensional and dimensionless parameters. 



Description 


Dimensional 


Dimensionless 


Hill radius 


Rk 


1 


Orbital frequency 


Q. 


1 


Headwind velocity 


"hv, 


4- 


Drag constant 


Kding 


KL, = 1/st 


Radius (proto)planet 


Rp 


ffp 


Stopping time 


t., 


St 


Collision rate 


Pcol 


p 


Impact radii" 


ha 


ba 



Note. — " For impact radii and velocities we intentionally use the same 
symbols, see also footnote|2l 



f sOoQ and the dimensionless headwind velocity, f„ 



velocity ^u, (Eq. (fTSll). the particle Stokes number St, and the 
protoplanet radius 

,1/3/ „ \-l/3 



Rp 

— =5.7x 10" 

Rh 



-3 1^ 



Ps 



3 g cm 



Iau/ 



(16) 



which mainly depends on semi-major axis aq- In Fig.|2]we give 
the relation between the dimensionless (u, and St to the physi- 
cal protoplanet size Rp and particle size s for two disk radii a. 
We have adopted Uhw = 30 m s"' and disk parameters that corre- 
spond (approximately) to a typical minimum mass solar nebula 
model (IWeidenschilling|ll977bl:lHavashi et al.lll985h . 

The full set of equations of motions in dimensionless form, 
dropping the primes, reads 



d_(v,, 

dt\vy 



3x/r^ 
2y.t - 3y/P 



2vy + 3x 



1 

Stlf^i 



■ C^ + 3.^/2 



(17) 



^u,^ 



ORh 



12.5 



gem 



-1/3 



/ Vh^ \l Rp Y ( a \ 
\30ms"i/\100km/ llAU/ 

(15) 
with Rp the radius of the (proto)planet. Note that due to the nor- 
malization to L)/,, (ui is primarily an indicator of the size of the 
(proto)planet rather than of the strength of the headwind ^hw as 
the latter is approximately constant throughout the disk. 



2.4. Dimensionless quantities 

Table [T] compiles some key quantities in both dimensional and 
dimensionless form. These include the dimensionless headwind 



The drag-free equations of motions are retrieved when St — > 
oo, which signifies that particles are not coupled to the gas. 
However, for Stokes number St < 1 particles are coupled to 
the gas and the im portance of the drag te rms becomes relevant 
or eve n dominant. iPetit & HenonI ( Il986l) and llda & Nakazawal 
(Il989h . working in the gas-free regime, only had to care about 
the first terms on the RHS of Eq. ( fTTI i and the equations of mo- 
tions did not include any parameter The addition of gas drag 
introduces two parameters: the velocity of the gas flow ^^, and 
the coupling parameter St. Together with the size of the planet. 
Op, these fully specify the problem; i.e., impact parameters b^- 
depend on these three dimensionless quantities only. Although 
not as clean as the drag-free equations, Eq. ( fTTI i still represent 
a significant reduction of the parameters involved (semi-major 
axis, particle size s, protoplanet size Rp, headwind velocity, gas 
density, etc.). In our parameter study we only have to care about 
these three parameters. 

We do not include the eccentricity in our prescription (and 
also not the inclination since the interaction is assumed to be 
2D). Rather, the initial velocity of the approaching particle is 
given by the radial drift equations for individual particles, ne- 
glecting the 2-body interaction term. These we will now review. 

3. The geometrical limiil 

Ignoring the 2-body interaction force, we will analytically solve 
for the particle's trajectory in the comoving frame. As we will 
soon see, impact parameters along a particle trajectory are gen- 
erally not conserved. We provide a general relation between the 
1/2 impact parameter at the interaction point (ba-) and its projected 
value on the x and y axes at any arbitrary point (Eq. (|2TI) ). This 
relation will be used later in Sect.|4|to obtain the impact radii b^. 



3.1. Steady-State velocities 

Without the two-body interaction term, t he motion of the particle 
fulfills the well-known drift equationsj Weidenschillindl 1 977al : 
iNakagawa et al..,1986^ Brauer et al. 200%: 



' Note that the Stokes number in this study simply indicates the di- 
mensionless friction time; it is not necessarily the same as the Stokes 



Vr = - 



f0 



2t;hwSt ^ 
1 + St- ' 

1 +St2' 



(18a) 
(18b) 



number used in turbulent studies, St = 
timescale of the largest eddies. For t^ 



dYoudin & Lithwickll2007h . 



arse 

m 



tjtt where tt is the turn-over 
= fi"' the definitions agree 



^ In this and the next two sections lengths (x,y,ba-, etc.) and 
velocities (v) are expressed in dimensionless (Hill) units, unless 
otherwise specified. 
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Fig. 3. Without the two body force, particle trajectories as wit- 
nessed from the comoving frame obey parabolas. Two trajecto- 
ries are shown: one that passes through the origin (yo{x)) and one 
that just hits the target. The corresponding impact parameter b is 
denoted by arrows. For curved trajectories b is not conserved due 
to the changing slope of the curves, here indicated by the angle 
0. At S we have that Ax^ - bs I sin 9s and Ai/s - bs / cos 9s ■ 
Ays is a conserved quantity. 



where v,- is the radial velocity and v^ the azimuthal velocity with 
respect to the local Keplerian rotation. Thus, in the context of a 
^xeii-rotating orbital frame we have v^ = y,. and Vy - v^- |Qo^, 
or in dimensionless units (divide by y;,) 



2^„St 



1+St' 



Cw 



l+St^ 



3 



(19a) 



(19b) 



and it can be verified that with these expressions the RHS of 
Eq. ( fTTl i vanishes when the two body interaction terms {-'ixjr' 
and -3y/r^) are omitted. 



3.2. The parabola solution 



Since 



dy_Vy^_J_ 3(1 + St^) 



(20) 



we immediately recognize that the particle's trajectory in the 
rotating frame obeys a parabola, y{x) - Ax^ + Bx + C with 
A = 3(1 -H St^)/8St^u„ B = l/2St, and C the integration con- 
stant, which is determined by the starting point S of the particle. 
We refer to the function that intersects the origin (C = 0) as 
yoix). The starting point ixos,ys) is defined to lie on this curve, 
ys = yoixos)^ see Fig. [3] where yoix) is plotted by the upper 
parabola. Another parabola with the same A and B (i.e., for the 
same particle properties St and ^u,) but with non-zero C is drawn 



in such a way that it just hits the 'target' at the origin. Its 'launch 
point' at y - ys is shifted over a length Ax^ . The vertical diff'er- 
ence between the curves, -C - Ays , is preserved. 

If the particle trajectories were straight, impact parameters 
would be the same everywhere in the (x, i/)-plane. However, due 
to the x--term this statement no longer holds for the general case 
of nonzero A. See Fig.[3j the impact parameter near the origin 
or at the interaction region, ba-, differs from that at S . Impact 
parameters are no longer conserved due to the change in dy/dx. 
The changing slope of the curve is indicated in Fig. [3] by the 
angle 9, where 9 is related to Eq. (l20b as = arctan(iij//iix). 

Using the properties of the parabola solution we can relate 
the quantities at / to those at 5 . At the interaction point / the 
impact parameter is b,r {=ap in the geometrical case) and the 
associated vertical width is Ayj - b^l cos 9i with 0/ the angle 
the parabola makes at this point. Similarly, Ays - bs I cos 9s 
and due to the invariance of Ay we therefore have that b^r - 
bs cos 0// cos 0s . 

The associated change in x at the starting point is Axs - 
Ays I tan 9s . This can be expressed in terms of the impact pa- 
rameter at the interaction point bg-, i.e.. 



Axs = 



^ys 



tan 9s cos 0/ tan 9s 



^b, 



V 



1 + {dy/dx)j 



(dy/dx) 



(21) 



where we used that 1 / cos 9 = 1 / cos[aTCtan{dy/dx)] = 
■\Jl + (dy/dx)-^. Of course, in the non-gravity limit we know al- 
ready that the impact parameter is just the planet radius, iio- = Up. 
However, Eq. (1211 1 is essential to interpret the numerical result of 
Sect. |4] That is, in our numerical integration (that includes the 2- 
body force) we will scan the x-axis for trajectories that lead to 
a collision and obtain a length scale Ax^ over which particles 
hit the target. Using the above equation we can then relate the 
obtained range of projected impact parameters at S , Axs , to the 
impact parameter at the interaction point, bg-. 

3.3. Collision rates 

The collision rate P in the 2D configuration is the product of 
the collision cross section, 2ba-, and the approach velocity, Ug = 

Jul + vl. At 5 we have that bsvs - Aysv^ - AxsVy. At / the 

collision rate equals bg-Va = apUg. Now, using Eq. ( |2TI ) we have 
that bg-U[ = bg-vj cos 9j = v^ tan 9s Axs = v^Ays and we see that 
the rates at / and S are equal and independent of the choice for 
the starting point ys. This result just reflects mass conservation. 
Thus, we obtain the collision rate: 



geo 



:4q'„ 



'(l+St^) 



(3V1+St')+4^J2 



1 + 



64St''f; 



■V2 



(22) 



where (dy/dx) has been evaluated at x = Up. The complexity 
of Eq. (l22T i may seem surprising for something as straightfor- 
ward as a geometrical sweepup. However, this is entirely due 
to the fact that Eq. ( l22l i covers several (velocity) regimes. In 
Appendix lAl we consider the asymptotic limits of Eq. (l22l i and 
show that these coiTespond to the expected sweepup rates (cross 
section x approach velocity) and to the findings of iKarv et al.l 
(Il993h . 
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Fig. 4. The minimum distance in units of Hill radii to the origin 
(center of the protoplanet), r^jn, as function of the x-coordinate 
of the starting point, xs (ys is fixed at 40 Hill radii). Plotted is 
'"min for a gas-free system {left) and a system that is characterized 
by the parameters St = 10 and ^^, = 1 {right). Several bands are 
labeled. The inclusion of gas drag shifts the bands to larger xs 
while merging several chaotic features within the chaotic c-band. 

4. Full 3-body integrations including gas drag and 
gravity 

4.1. Description of tlie adopted algoritlim 

We perform a parameter study of Eq. ([TtT i. varying fu, and St. 
For the dimensionless headwind velocity runs were performed 
at (^ = 0.01, 0.03, 0.1, ... 10'' and for the Stokes number values 
of St = lO-'', 3 X 10-* . . . lO'' were sampled. Thus, we obtain a 
grid ofl7x 13 - 221 different combination of ^u, and St. Not 
every combination is equally likely. Indeed, the parameter space 
samples areas where our key approximations (linear drag law, 
constant gas density) lose validity, but we intentionally sample 
a broad range of values to verify the validity of our analytical 
expressions (see Sect.|5ll. 

For each combination of ^^, and St we numerically deter- 
mine the function rmin{xs). Particles are launched from a start- 
ing point {xs,ys) where ys is fixed and xs is varied, see Fig.|3] 
The initial velocities are given by Eq. ( fT9] l. Depending on the 
sign of v,j{xs) the initial (/-position (+ or -) is determined, such 
that the initial motion in y is always directed tow ards the planet. 
Here, we fix \ys\ at 40 (llda & Nakazawal [19891) . Particles that 
leave the computational domain (when \y\ > 40 or x < -40) are 
no longer followed. For a certain xs we numerically integrate 
Eq. ( fTTl i adopting an relative error of at most 10-^. As our in- 
tegrator we use a fifth-order Runge-Kutta sche me with timestep 
control (lFehlberg|ll969t IShampine et al.lll979b . After the calcu- 
lation has terminated we determine (and store) the minimum dis- 
tance, rniin. In this way r„im{xs) is obtained, see Fig.|4] Projected 
impact parameters Axs are then obtained from the r^in{xs ) curve 
by summation over the intervals where rmin < ap, i.e., 

Axs = I dxsH{ap-r^in[xs]), 
where H{t) is the Heaviside step function. 



(23) 



H{t) 



1 r>0; 
f <0. 



(24) 



(Simply put: we only include the orbits that hit the target). 



As r^in is occasionally found to vary steeply with xs, fine 
sampling of the x-axis is required. Therefore, we sample our pa- 
rameters space {xs) adaptively. We start out with intervals of 1 
Hill radii, e.g., {xs} - 0, 1, 2, 3, In the next level the inter- 
val spacing is reduced by a factor 10, S = 0.1. However, we 
only treat the points that fulfill the condition r^in{xs) < F6, 
where F is empirically fixed at lO"* (see below). For example, 
if '"min(O) - 300 this point will be skipped in the next iteration 
of the algorithm. If xs fulfills the condition, however, then both 
sides will be scanned; e.g., if rn,in(2) = 15 < Q.IF then, we will 
additionally perform calculations for xs - 1.1, 1.2, ... 1.9 and 
Xs = 2.1, 2.2, . . .2.9. In this way we reduce the number of cal- 
culations but are still able to obtain a good assessment of Axs 
for low a p. 

Despite this optimization, we were forced to perform a rela- 
tively large number of integrations, i.e., a large F. The reason is 
the presence of very narrow, chaotic bands. In Fig. |4]band b near 
Xs = 2.0 is a regular band since rmin varies smoothly with xs. 
A low F value suffices to pick up this feature. However, bands 
a and d are very narrow and show (if one would zoom in) addi- 
tional substructure. These chaotic bands are not resolved when 
choosing a low F. In fact, there is no guarantee that our algo- 
rithm will pick up every band since they can be very narrow. 
However, with F - 10^ we do obtain a good correspondenc e 
to previous works ('Pet it & Henorf l986Vld a & Nakazawal 19891) . 
also matching the substructure within the narrow bands shown in 
Fig. HI 

Following the discussion in Sect. [3] we emphasize again that 
the starting points {xs values) are not fundamental, but depend 
on the choice of the starting point ys . Taking a different value 
of ys, e.g., ys = 80, will shift the features of Fig.|4j5 towards 
higher xs values. In addition, the spacing (width of the features) 
will be different. The only conserved (physical) quantity is the 
mass flux, i.e., the integral of J v,j{xs)dxs over the width of the 
feature - independent of the choice of ys . 



4.2. Orbits including gas drag 

Figure |5] provides several examples of particle trajectories that 
experience gas drag. Figure|5h shows the trajectories for f^, = 1 
and St = 10 with different starting points xs (Figure |4] contains 
the same parameters). For a 'standard' nebula setting, these pa- 
rameters correspond to ~m-size particles accreting onto a alO-' 
km planet, see Fig.|2] Due to the large Stokes number, the in- 
fluence of the gas is relatively weak and the orbits bear a close 
resemblance to the gas-free, three body regime. The xs =3.8 
trajectory enters the Hill sphere, where it experiences a close en- 
counter, at Tmin = 2.8 X 10-2, before leaving the Hill sphere. It 
then re-emerges later at negative x due to the combined effects 
of inwards radial drift and Keplerian shear However, the particle 
that started out at xs =3.9 is captured within the Hill sphere and 
experiences strong orbital decay due to gas drag. 

Figure |5}' shows orbits for smaller particles of Stokes num- 
ber 0.01. Four orbits are shown of which two lead to accretion. 
Clearly, there is a contest between the gravitational pull of the 
planet and the aerodynamic pull of the gas flow. Once close 
enough, gravity always wins. All orbits with 0.38 < xs < 0.74 
are accreted; there are no close encounters. Accretion is inde- 
pendent of the physical proportion or internal density of the pro- 
toplanet; once a particle's angular velocity about the planet is 
damped by drag, it settles radially at its terminal velocity. The 
only relevant physical quantity is the mass. This mode of accre- 
tion reflects the capture mechanism of Fig.|5t. We will refer to 
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Fig. 5. Examples of planet-particle interactions for different values of the dimensionless headwind velocity (u> and coupling pa- 
rameter St. For typical nebula parameters particles of St = 10 correspond to loosely coupled m-size particles, whereas St = 0.01 
are more strongly coupled cm-size particles, see Fig. [2] Likewise, f„, = 1 corresponds to protoplanets of Rp ~ 10^ km in radius, 
while ^u, = 100 corresponds to Rp ~ 10 km planetesimals. (A) Two particles of St = 10 experience a close encounter within the 
Hill sphere (dotted circle). The xs - 3.9 particle is captures and settles to the planet, whereas the other particle is ejected from 
the Hill sphere (The Keplerian shear eventually causes it to resurface at the other side of the Hill sphere). (B) Strong gas coupling, 
St = 0.01. There is a competition between the gravitational pull of the planet and the drag force directed towards negative y. (C) 
Close encounters at large f^ without settling (see inset). (D,E) Examples of particle trajectories originating from interior orbits. (F) 
Radially approaching orbits. 



orbits like the xs = 3.9 curve in Fig.|5^ as gas drag induced or- 
bital decay, whereas the accretion mode in Fig.jSj" is referred to 
as settling and draw the dividing line at St = 1 . 

On the other hand Fig. [5]:, which features a larger dimension- 
less headwind (meaning: a smaller protoplanet) of ^u, = 100, 
does not display the settling behavior. Here, particles can only 
be accreted due to the finite size of the target. The xs - 0.796 
trajectory has a minimum distance of rmin = 5.0 x 10 "*; the 
Xs = 0.8 trajectory r^in - 4.5 x 10""'. Clearly, for a planet size 
Op ^ 1 the impact parameter in Fig.|5}; is much less than for 
the settling orbits of Fig.|5]5. Since the Stokes numbers are the 
same, the reason must be due to the larger headwind velocity f„,. 
This is understandable since particles of St <sc 1 approach at the 
headwind velocity (ua x ^^) and a large Vg is not conducive for 
accretion. 

In the lower panels of Fig. |5] we vary either the Stokes num- 
ber (particle size) or ^^, (protoplanet size) with respect to the 
panel above. For a Stokes number of lO"*, see Fig.|5}l, the ef- 
fects of gas-drag are even less pronounced and it becomes more 
difficult to capture these (big) particles within the Hill sphere. 
Moreover, if such a particle would be captured, it takes longer to 
finally accrete this particle due to orbital decay. Another differ- 



ence with Fig. |5^ is that the St - lO"* particles can now also enter 
the Hill sphere from interior orbits (negative xs). In Fig.|5^ the 
strong radial drift still prevents particles from entering the Hill 
sphere from the negative {/-direction; however, for St = 10^ the 
radial drift is sufficiently reduced to render the situation more 
akin to the symmetric gas-free limit. 

The ^ui - 0.1 orbits in Fig.|5^ also feature accretion from 
particles approaching the planet from interior orbits, which the 
^tt, = 1.0 orbits of Fig.|5j' were not capable of. The dimension- 
less headwind parameter of ^„ = 0.1 corresponds to a very 
big planet (in the canonical model) for which, as we will dis- 
cuss below, the constant gas density background is unrealis- 
tic. Alternatively, it can represent a smaller protoplanet in a 
nebula where the dimensional headwind is, for some reason, 
strongly reduced. In any case, we see that low ^^ tends to make 
the interactions more symmetric. This can be seen from the 
Eq. (I19bl i: low f^, or large St reduce the contribution from the 
non-symmetric headwind term, ^u^/(l + St^). 

Figure |5f shows, however, that for St = 10 and (^, = 100 
the picture is anything but symmetric. The particles approach 
the planet from a very radial direction (jc-direction) - at least, as 
seen from the perspective of the planet. The point is here that 
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both St and ^(„ are large. Thus, both planet and particle move 
at a Keplerian velocity (in the azimuthal direction) but, due to 
the large (,^„ the particle still suffers a significant radial drift, 
which outweighs the effects of the Keplerian shear As a result, 
the situation is similar to Fig.|5];: accretion does only proceed 
through close encounters. 

4.3. Collision rates 

We obtain the (dimensionless) collision rate from the encounters 
that hit the protoplanet, i.e.. 



J 



P{a,„SU(u,)^2 dxs\uJxs)\H(ap-r^in[xs]), 



(25) 



with Vy{x) given by Eq. (I19bb and H{t) the Heaviside step func- 
tion. Figure |6] plots contours of P(^u„St) for a planet size of 
ap = 10"^, which corresponds to an heliocentric distance of ^5 
AU. The reader must realize that P is expressed in dimensionless 
units; the large rates that can be seen at large (ui (and St < 1) are 
less impressive upon multiplication by R^Vh k ^~^ (see Eq. (fTSl)). 
In fact, these high P values are consistent with the geometrical 
sweepup rates of Eq. ( l22l i. However, the expression in terms of 
dimensionless units is useful since we can directly compare it to 
the gas-free limit for which Pgf ^ llaj" — 0.35, see Sect. 12. II 
For large St and small (,^, P converges to Pgf, the expected be- 
havior. However, for the remainder P deviates significantly from 
Pof. We sum up the main features: 

1. Particles of St ~ 1 accrete very well when the headwind 
velocity is low. There is a distinct peak at (St, ^u,) = (0, 0); the 
accretion rate is here 20 times larger than Eq. ([TJ. However, 
there is a very sharp transition between 1 < ^^ < 10. 

2. For large ^u,, P - P„^o is larger than Pgf although no gravita- 
tional focusing takes place. The sweepup is so effective due 
to the strong headwind. 

3. For St > 1, the band St ~ ^^ features a maximum in P. 

4. For St «; 1 and low ^u, (large planets), collision rates are 
lower than Pgf. Tiny dust particles stay 'glued' to the gas 
due to their strong coupling, preventing accretion. 



5. A simple model for gravo-gas interactions 

5.1. Model outline 

We present a simple model for the impact parameter ba-. The 
model is summarized in Fig.|7l where the three relevant regimes 
for the impact parameter are shown. In the hyperbolic regime 
encounters follow the two-body approximation. The usual grav- 
itational focusing formula applies. Keplerian shear is unimpor- 
tant. In the settling regime particles settle to the target and the 
impact parameter is independent of the planet size, a p. For this 
reason, impact parameters can become rather large. Finally, in 
the three-body regime, the encounter proceeds along the lines 
of the drag-free three body encounters at low energy. However, 
the presence of the gas now causes some particles to be captured 
within the Hill sphere; these orbits decay and this enhances the 
accretion rate. 

5.1 .1 . Importance of three body encounters. 

It is clear that for St » 1 the encounter cannot be described 
by a 2-body interaction, but should include the solar gravity. 
But what is the transition between the 2-body and the 3-body 
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Fig. 7. Illustration of the three accretion regimes. In the hyper- 
bolic regime interactions are 2-body encounters and the standard 
formula for gravitational focusing applies. In the settling regime, 
accretion proceeds through settling, which enhances the impact 
parameter bg-. Else, for St > max(^u,, 1) the solar gravity affects 
the encounter and b^r is increased with respect to the gas-free 
limit {b„f, Eq. (|5]l) due to particle capture within the Hill sphere. 
The gray band approximately indicates the zone where the set- 
tling solutions permeate into the hyperbolic regime. The letters 
a-f correspond to the parameters in the panels of Fig. |5] 



regime in the presence of gas drag? A passage through the Hill 
sphere typically takes a time of the order of the orbital period. 
Thus, at first sight, we can draw the boundary at St = 1 since 
particles of lower stopping time will be strongly affected by the 
gas. However, a large headwind velocity (ui has the same effect. 
Particles that experience a drag force ^m/St > 1 are blown out of 
the Hill sphere. Thus, 3-body effects are reduced to the region of 
parameter space where St > 1 and St > (ui, see Fig.|7] 
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5.1.2. Two body regime: settling- and liyperbolic interactions 

We consider an interaction at impact parameter b - ba- at an 
approach velocity u^. The strength of the gravitational force is 
fg = 3/b~ and the interaction timescale, f„ ^ b/ua- The lat- 
ter quantity can be compared to the particle's response time 
St. When ta < St gas drag can be neglected during the en- 
counter; the induced velocity change is Au = fgta - 3b/ua. 
However, if ta > St the particle's velocity equilibrates towards 
Av = A^set = fgSt = 3St/b-. 

For St < 1 the approach velocity can be approximately writ- 
ten as Va - 3b/2 + (u,. For low f^ we therefore can expect settling 
since encounter timescales ?„ are long. However, for large ^u, set- 
tling will be prohibited: either the induced change Av is too little 
(at large b) or the interaction timescale too short for the particle 
to obtain its settling velocity (at low b). 

To see this quantitatively, the minimum impact parameter for 
settling is b* - ^^jSt and the corresponding velocity change is 
Ay* - 3/^^St. In order for the particle to settle to the central ob- 
ject, the direction of the particle has to change over a large angle, 
i.e., Av ~ Va- In fact, we obtain a better correspondence with our 
numerical result if we set the required velocity change to ValA-. 
Analytically, it can be shown that this is the required change for 
St <s; 1, see AppendixlB] Thus, accretion through settling takes 
place when Val^ < Anset ^ Ay* and disappears when y„/4 > Ay*. 
At the boundary between the settling and hyperbolic regime it 
is allowed to take Vg - (,w (as can be verified a posteriori). We 
then have that for f^St > 12 settling is no longer possible, corre- 
sponding to a critical Stokes number 

St* = i^, (26) 

above which accretion through settling will no longer occur. 



5.1.3. Tine settling regime 

Assuming the settling regime, particles at impact parameter b 
experience a velocity impulse of Ayset = 3St/fe^, which should 
equal ValA- for accretion. Since y,, ^ 3bl2+(ui in this regime (note 
that we do not neglect the shear term since it becomes important 
at low (ui) the condition Ayjet = ValA- require us to solve the cubic 
equation 



b' + ^b^ ■ 



8St = 0. 



(27) 



The (real, positive) solution to this equation is denoted b^, 



5.1.4. Hyperbolic regime 

At large f^, the encounter is fast and the presence of gas drag can 
be ignored during the encounter. For accretion we now require 
(by conservation of angular momentum) that Ay = 3b/ua > yesc 
which is much larger than in the settling regime. For the im- 
pact radius we can just take the standard expression of the 
gravitationally-enhanced cross section. 



'hyp 



1 + 



^esc 




(28) 



where ygsc = ^j^|ap in Hill units. For the approach velocity v„ 
we now also include the horizontal velocity components (i.e., 
Vx) since for St > 1 y^ becomes dominant. However, it is fine to 



neglect the shear term in Eq. ( I19ab since /^hyp "^ 1 ■ Thus, for the 
approach velocity in the hyperbolic regime we can write 



'^ajhyp ~ ^u 



5.1 .5. The three body regime 



Vl + 4St2 

l-HSt^ 



(29) 



Without gas drag the effective impact parameter for collisions 
is Eq. (|5]l, /jgf = l.Vttp , in dimensionless units (see Sect. 12. lb . 
Gas drag adds another component to the impact parameter on top 
of Eq. (|5]). Figure |4] neatly illustrates this behavior. The chaotic 
band c in the gas-drag simulation has collapsed. Particles en- 
tering at the corresponding xs -values are captured and decay to 
the central object on a timescale ~St. If the gas inside the Hill 
sphere is removed within this timescale, these particles will be- 
come satellites; however, here we will simply assume that all 
captured particles contribute to the collision rate. 

Because the accretion in the dissipative 3-body regime is 
determined by the behavior of the chaotic zones, it is difficult 
to provide an analytic model for the enhanced b^. The chaotic 
zones are especially susceptible to collapse, because these par- 
ticle trajectories are characterized by many revolutions, trough 
which a lot of energy can be dissipated. In the gas-free situation 
one requires a positive energy J to enter the Hill sphere. 



,1,33.9 

J — —V x" H — ; 

2 r 2 2 



(30) 



and once J becomes negative in the Hi ll sphere, e.g., b y inelastic 
collisions, the body becomes trapped (IOhtsukil|l993h . Here, we 
face a similar situation where the gas drag is responsible for the 
energy removal. Unfortunately, in our case an analysis in terms 
of the Jacobian is not so meaningful as the gas flow can also add 
energy; i.e., J is not conserved and bodies with 7 < can still be 
'blown out' of the Hill sphere. However, the picture - that gas 
drag can trap particles - is still the key. 

Empirically, we find that the impact radius is increased by 
a term 1 /St, which corresponds to the dissipated energy over a 
revolution. For these reasons, we add a term proportional to 1/St 
to Eq. ©, 

1.0 ,n 1-0 

^3b = %+— = l-7f^y +— ' (31) 



St 



St' 



where the 1 .0 constant is obtained empirically. As we will see in 
the next section Eq. ( ISTT i fits the general trend well, but it cannot 
reproduce the impact radius at every Stokes value. 

5.2. Comparison to numerical results and fine tuning of the 
recipe 

Figure |8]compares the impact radii obtained from the numerical 
integrations (symbols) with the analytical prescriptions (curves) 
for a headwind velocity of (^i = 10 and for a planet size of 
Op - lO""* and 10"''. We have used Eq. (|2TI ) to convert the pro- 
jected impact parameters Axs to true impact parameters at the 
interaction point (ba-). To do so we used the parabola solution, 
Eq. (l20t . to evaluate the gradients (dy/dx) at the starting point 
(xs) and at the interaction point xj. We determine the maxi- 
mum value of Xs that resulted in a collision with the planet at 
the specified Up and took this value to compute {dy/dx)s. To 
compute {dy/dx)i we evaluated Eq. (l20t at the approach radius 
^app. Here, for ^app we took the impact radius obtained from our 
analytical model described above, except for interactions in the 
3-body regime where we always use /japp = 2.5. 
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Fig. 8. Impact radii from the numerical integrations (symbols) 
and analytic fits (curves) for a headwind velocity of ^„, = 10. 
Vertical lines St* and ^u, distinguish the settling, hyperbolic, 
and three-body regimes. Analytic fits from our recipe are shown 
by the dashed-black, solid-gray, and solid black curves, respec- 
tively, and denoted Sset, ^hyp and i>3b. The corresponding im- 
pact radii obtained from the numerical integrations are shown 
by crosses (Op - lO""') and circles (ap = 10"^). 



At low St the interactions take place in the settling regime. 
Impact radii are rather large, particularly near the St* = 12/^^, 
transition line, and independent of Up (the cross and circle sym- 
bols overlap), implying that the physical impact parameter is 
larger at larger disk radii. For intermediate Stokes numbers the 
hyperbolic regime is valid and impact radii are much smaller. 
However, for St > ^^, impact radii once again increase. The 
behavior is rather erratic, though, with peaks at St - 10 and 
300 and a depression at 10^, valid for both ap. We found that 
this complex behavior can be attributed to the trajectories that 
originate from the third quadrant. Initially, for low Stokes num- 
bers, these are absent due to the strong radial drift. However, at 
a critical Stokes number the contribution of particles approach- 
ing from interior orbits (negative ys) becomes important. We do 
not have a full understanding how these outliers can be modeled 
analytically. 

The analytic fits to the various regimes are given by the 
dashed curve (for settling), solid gray curve (hyperbolic) and 
solid black curve (three-body). From Fig.[8]it is obvious that the 
transition between the settling and the hyperbolic regime is not 
so sharp. Even particles that have St > St* display settling be- 
havior. For these reasons, we have extended the validity of the 
settling regime beyond St* by adding an exponential term, i.e.. 



^,se, = ^seteXp[-(St/St*)n, 



(32) 



where b^st is the solution of Eq. dZTJ i and y a constant that we 
empirically fix at y = 0.65. The impact parameter in the hyper- 
bolic regime, St* < St < C,w, is then given by the maximum of 
^set and /jhyp- The little depression that can be seen at St = 1 is 
caused by the fact that the approach velocity Eq. (|29] l has a maxi- 
mum here. For larger St the approach velocity strongly decreases 
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Fig. 9, Contours of log P according to the analytic prescription. 
Curves are the same as in Fig. |6] 



and the gravitational focusing strongly increases. Nevertheless, 
impact parameters at St = (^, are even larger than ^hyp and are 
better fitted by b^.^. Initially, gas drag very effectively captures 
bodies within the Hill sphere and impact radii are rather large. 
However, the capture probability decreases as 1/St and for large 
St we retrieve the gas-free limit, Eq. (|5]l. 



5.3. Collision rates 

In Fig.|9]we plot contours of the collision rate, that is, we plot 
P = 2bcrVa as function of St and ^„, using the prescription out- 
lined in Table |2] This figure should be compared with Fig.|6] 

The curves in Fig. |9] are much smoother due to the much 
finer grid that the analytic formulation permits. The transition 
Unes, St - St* and St = (w are clearly identified. Our analytic 
formalism fails to reproduce the log P - -0.5 band towards the 
upper-right of Fig.|6l However, the overall match is satisfactory; 
for 90% of the 22 1 grid points the analytic and numerical results 
lie within 30% of each other 



5.4. Summary of impact parameter recipe 

Table|2] provides an executive summary of how the collisional 
parameters can be obtained using the analytic prescription. 
First, one converts the physical parameters (headwind velocity, 
disk radius, friction time, etc.) into the dimensionless quanti- 
ties ^u], o:p and St. The corresponding impact radii for the three 
regimes are calculated in the second step. Then, in step 3, the 
appropriate collision regime is determined by comparing the 
Stokes number with St* (Eq. ( |26] |) and ^^ (Eq. (flSll). see also 
Fig- ID Dependent on the applicable regime, the final impact pa- 
rameter is obtained by taking the maximum of two impact radii 
(step 4). Other quantities {ua and ^app) also depend on the colli- 
sion regime. 

Then, these results can be converted back to physical units 
by multiplication of Ri, and vi, = Ri,Q. (Eq. ^) for, respectively, 
lengths and velocities. The 2D-collision rate is then obtained 
from Eq. (|3]l. The 3D-collision rate may be estimated by mul- 
tiplication by a factor max{l, Hp/bg-) (see Sect.|6]l, where Hp is 
scaleheight of the particles. 
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1. Calculate dimensionless parameters: 


(ui (headwind velocity) 

Up (planet size) 

St = fjQ (Stokes number) 


Eq.dBl 
Eq.(fT5l 
Eq.ll2j 




2. Calculate impact radii: 


bhyp 

bib 


Eq. m, Eq. (03 

Eq.m 
Eq.(|3l) 




3. Determine Regime: 


St<min(l,12/<r^,) 
Settling 


Hyperbolic 


St > max(f„„ 1) 
Three body 


4. Results 

Impact radius (accretion), ba-: 
Approach velocity j)„: 
Approach radius b-^pp-. 


max(^sct,Vo) 

36^/2 + 4, 
bo- 


max(b,a, bhyp) 
Eq. ^29J 
b^ 


max(i3b,fegco) 

3.2 
2.5 



Note. — Description of impact radii: b^^o, geometrical impact radius (= or,,); b^st impact radius in settling regime; b^ct, modified b^^i (to cover the 
transition regime); fchyp impact radius in the hyperbolic regime; bjb drag-enhanced impact radius for the 3-body regimes; b^pp, approach distance. 



6. Significance to the growth of pre-planetary 
bodies 

In the previous sections we have outlined a general approach to 
analytically derive impact radii and collision rates in Hill coor- 
dinates. But what does all of this imply for the growth of pre- 
planetary bodies? Perhaps the best way to illustrate this point is 
to calculate the accretion timescale 



j2D _ 



M ^np.Rlli 



dM/dt fcoiS 



32P 



D. 



(33) 



where P^oi is the dimensional accretion rate and P - PcoilRh^h 
the dimensionless, the quantity plotted in Figs. |6]and|9l From 
Eq. ( |33] | we see that the accretion timescale is inversely propor- 
tional to P but also scales with Rp. We further write Eq. ( l33l l in 
terms of semi-major axis oq by substitution of Eq. ( fTSI l for Op 
and Q.{a) for a solar-mass star 



n2D 



J.U ^ 6^ / Ps 

P \3gcm- 



-1/3 



1 g cm 



Rn 



100 km/ 1 AU/ 



-1/2 



yr. 



(34) 

The inverse dependence on disk radius may seem surprising but 
one has to realize that P via Op and (u) implicitly depends on 
a. Nevertheless, Eq. ( l34l i shows that 2D accretion may be espe- 
cially advantageous in the outer disks. 

The 2D regime, however, may not be applicable to small par- 
ticles since any breath of turbulence will stir them up. The height 
of the particle layer may be obtained by equating particle diffu- 
sion and settling timescale; i.e.. 



(35) 



— X min 1 , 1 / — 

H, r V St 



dDubrulle et al.l 119951 : ICarbalUdo et"ai1 120061: 

I Youdin & Lithwickl l2007h where a, is the IShakura & Sunvaevl 
( Il973h viscosity parameter for turbulent diffusion. In a 3D 
setting, the particle scaleheight can exceed the impact parameter 
by, then, only a fraction, b,r/Hp, of the particles take part in 
the interaction and the accretion timescale is correspondingly 
longer. 



j-'^U ^ yZ 



X max 



1,^ 



(36) 



(with bcr in physical units). The 3D correction factor signifi- 
cantly increases collision timescales for small particles (Hp is 
large) and the hyperbolic regime (ba- is small). 



In Fig.fTOh we have plotted contours of the 3D growth 
timescales for a disk radius aq = 5.2 AU, p, = 3 g cm"-* (making 
Op = 10-^^), 2 = 2 g cm-2, Hg = 0.25 AU, A^hw = 30 m s"', 
and a, - 10 ■*. The calculated accretion timescales assume that 
all the solid density is in particles of a single size. Due the inclu- 
sion of the Hp/bij- factor the structure is quite different from that 
of Fig.|9l However, the contrast between the hyperbolic and set- 
tling regimes is still clearly visible and has in fact even increased 
due to the correction factor for the vertical structure. Note that 
for the bigger bodies, which settle into a thin plane, Fig.fTOlstill 
assumes that their eccentricities and inclinations are absent (low 
velocity regime). 

Figure[TO] tells a few interesting points. First, it can be 
clearly seen that growth of ~km-size planetesimals by accre- 
tion of small particles (St < 1) takes a (perhaps prohibitively) 
long time. Two mechanisms conspire. First, the small parti- 
cles couple effectively to the gas which dilutes their number 
densities near the midplane where the planetesimals are resid- 
ing. Of course, this statement depends on the strength of the 
turbulence that prevents the particles from settling effectively: 
timescales will be shorter for lower turbulent strength parame- 
ter, a,. However, even in a compl etely laminar disk we m ay ex- 
pect shear turbulence to develop (IWeidenschillind | l980|). which 
strength may be equivalent to arvalue s of ~ 1 0"^ (iJohansen et al.l 
120061: ICuzzi & WeidenschiUindl2(306h . Second, small particles, 
being strongly coupled, move with the gas, at a relative velocity 
of Au X yjiw s> i)esc, where (;},„ is the velocity of the headwind 
and Liesc the escape velocity of the planetesimal. Therefore, small 
particles lack gravitational focusing and it is hard to avoid the 
conclusion that sweepup of small particles by ~km size plan- 
etesimals is slow. In order to grow, planetesimals have to accrete 
among themselves. 

However, the situation completely reverses when protoplan- 
ets sizes reach ~10^ km: for these bodies, accretion of cm to 
m-size particles becomes very rapid: in only ~10^ yr the proto- 
planet can double in size. This is entirely due to the increased 
cross section in the settling regime. For 'optimal' parameters 
(St ~ I, im ^ 1) the combined effect of gravitational focusing 
and gas damping results in impact parameters of ba- ~ Q.SRh 
- larger than what hitherto has been thought possible (fegf, see 
Eq. (|5]l). Accreting at impact parameters of the order of the Hill 
sphere is fast in any case but since Ri, increases with disk radii it 
is especially impressive for the outer disk, see Fig.fTOb. 

For St ~ 1 particles, accretion is fast - even though it is 
inefficient due to the strong radial drift. We denote the probabil- 
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Fig. 10. The 3D growth timescale T^^ as function of Stokes number (or particle size) and size of the protoplanet. Contour lines of 
logjQ T^^/yi are shown. All the solid density is assumed to be contained in particles of the indicated size. (A) At 5 AU for a surface 
density of 2 = 2 g cm"^. (B) At 35 AU for a surface density of Z = 0.1 g cm"^. In (A) the wider spacing between the tickmarks of 
the upper axis indicates particles enter the Stokes drag regime. 



ity that radially-inward drifting particles become accreted by the 7. DisCUSSion 



protoplanet Peft- Since the drift flow is Inavrll (a:0.15Me yr 
for St = 1 particles at 5 AU) Peff is given as 



eff 



fcol 

2nav,- 



X max 1 , — 



(37) 



(where we have again included the scaleheight correction fac- 
tor). For a protoplanet of Rp - 10^ km we find that St ~ 1 
particles are accreted at an efficiency of only PefF ~ 10""*. The 
efficiency increases away from the St = 1 line [u,- is lower) and 
towards larger protoplanet size and larger disk radii (larger R/,). 
[Levison et al. (2010), using N-body techniques, also finds that 
accretion of small fragments was very inefficient (and concluded 
that it was therefore unlikely) except for a few specific particle 
sizes, that may have coincided with the peaks in Fig.fTO] To in- 
crease the accretion efficiency, smaller particles have to be ac- 
creted. Dust fragmen tation ( Birnstiel et al. 200a ) or radial dif- 
fusion dCieslal l2009h may be promising mechanisms to retain 
mm-size particles in th e (outer) disk, where t hey are observed on 
~Myr timescales ('e.g.. lLommen et aDl2009l : iRicci et alJl2010h . 

It is instructive to compare the accretion timescales of 
Fig.fTOk to detailed hydrodynamical simulations involvingSt ~ 
1 particles ( Johansen et al. 2007; Johansen & Lacerda 20 fy)- In 
Ijohansen et al.l (l2007h a dense particle layer of St ~ 1 boulders 
collapses into a Ceres-mass planet (R ~ 500 km), that rapidly 
accretes the remaining boulders on timescales of perhaps 10 yr 
Although from Fig.fTOh a Ceres-mass protoplanet in combina- 
tions with St ~ 1 particles form the optimal growth conditions, 
our accretion timescale of 10^ yr is still t wo orders of magnitud e 
higher than what can be inferred from iJohansen et al.l (l2007h . 
Howe ver, a direct comparison is perhaps not so meaningful since 
in the iJohansen et al.l (l2007h simulations the St ~ 1 particles 
are highly clumped and exert a strong feedback effect on the 
gas (IJohansen & Lacerdall20 Id discuss some alternate settings). 
Feedback effects are not taken into account in this study. 



7.1. Neglected effects 

We summarize the key assumptions that have been employed in 
this study: 

- a drag law linear in velocity; 

- neglect of resonance trapping of particles; 

- a smooth, laminar disk (only drift motions) without local 
pressure fluctuations; 

- neglect of the gas flow around the protoplanet and of a pos- 
sible atmosphere surrounding the proto(planet); 

- a dynamically cold protoplanet on a circular, non-migrating, 
orbit. 

The assumption of a linear drag law implies that this study - and 
in particular the analytical prescriptions that have been derived 
- apply for particles smaller than ~iniax (see Eq. ( fTlT l) only. But 
this still covers an appreciable size range, especially for the outer 
disk. 

The adopted flow pattern in our study is unrealistic since it 
does not take account of the presence of the planet. Of course, 
streamlines will have to bend around the object and this will af- 
fect the motion of the particle. Tiny dust grains can only col- 
lide with a dust aggregate when the aggreg ate size is less than 
the mean free path of the gas molecules dWurm et al.l 120011 ; 
ISekiva & Takedall2005h . However, this restriction probably ap- 
plies only for small particles. If we assume that the flowlines 
change over the lengthscale of the protoplanet, it follows that 
particles of fj < Rplv^^ are too tightly coupled to the gas to be- 
come accreted. For example, for Rp = 100 km only St < 10"^ 
particles are affected and our model overestimates the (already 
low) accretion rates. More serious, perhaps, is our neglect of col- 
lective effects due to strong particle volume densities, as this will 
provide a feedback effect on the gas, affecting both the flow pat- 
tern as well as the drift rates. This will probably be important for 
St ~ 1 particles and could significantly enhance the accretion 
rate (see our discussion at the end of Sect. |6]l. 

For simplicity, our analysis only included drift motions. We 
have, for example, neglected a systematic accretion (or decre- 
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tion) flow of the gas. In the turbulent a-model this gas moves in 
at a velocity ~a,CgHg/a. Equating this expression with the ra- 
dial drift velocity (Eq. ( |19at ) we find that for particles St < a, 
the systematic accretion flow will dominate. This will affect the 
expressions for the accretion rates. Likewise, turbulent motions 
can become more effective to move particles around than drift 
motions, which aff'ects the input parameter ;;„ in our model de- 
scribed in Sect. |5] In the ff,-model large eddies transport parti- 
cles at velocities ~a, Cg (ICuzzi & Weidenschillindf2006l) and 
turbulent velocities will dominate the drift motions for a, > rj 
(see Eq. (fTOt). The presented model may still be valid though, 
if the turbulent motions are included in the definition of the ap- 
proach velocity, Va- 

Mean motion reson ances may halt the particle long before it 
drifts to the Hill sphere ( IWeidenschilling & Davisll985h . In such 
a situation the inward-directed drag force is balanced by t he out- 
ward resonant perturbations. IWeidenschilling & DavisI (Il985h 
showed that the strength of the perturbations is proportional to 
the planet's mass and to the resonance number j. Thus, smaller 
particles (which experience a stronger drag force) move into a 
higher resonance. However, this trend will not pursue indefi- 
nitely as at s ome maximum / reso nances will overlap and the 
effect is lost. iPaardekoopej (l2007h simulated particle accretion 
onto gas giants and showed that ~m-size particles (St ~ 1) avoid 
resonance trapping for Jupiter-mass (Mj) planet. For a planet of 
0.1 Mj the critical Stokes number has risen to St = 10 and for 
lower mass planets it will even be larger. Therefore, our results 
are not so much affected by resonance trapping when the proto- 
planet and core-formation stages are considered. 

Our assumptions of a completely 'inert' protoplanet is also 
peculiar In the oligarchic growth regime (where most of the 
mass resides in leftover planetesimals) dynamical friction will 
keep the eccentriciti es of the most massive bodies small (e.g., 
iKokubo & Idal l2000h : however, if most of the mass is trans- 
ferre d to the protoplanets t heir motion will become eccentric 
(see iKarv & Lissaueij 1 19951 for accretion probabilities of pro- 
toplanets on eccent ric orbits). Likewise, type-I radial migration 
dTanaka et alj2002h is not incorporated in our framework . These 
effects will again become important for already evolved proto- 
planets of masses > O.IM®. 

The feedback effect of the protoplanet on the structure of 
the gas disk is also neglected. The protoplanet's gravity influ- 
ences the gas disk at larger distances, which could invalidate 
our approximation of a smooth (global) pressure. Indeed, par- 
ticles have a tendency to drift to high pressures regions and the 
St ~ 1 particles may be most affected by this process, piling 
up at a pressure bum p instead of proceeding to the protoplanet. 
iPaardekoopeJ (l2007h finds that this effect (together with the res- 
onant trapping of bigger bodies) shuts off all accretion of parti- 
cles sizes above ~10 - 100 /^m! Howe ver, this particle trapping 
is applicable for evolved planets only. iMuto & Inutsukal (l2009l) 
found that the criterion for particle trapping 



H, 



Mp > i]^M^ ~ 10 Me 



(38) 



for solar mass stars (Equation ( |38] ) is apparently independent of 
particle size or Stokes number). 

Long before this size is reached, protoplanets bind the 
nebular gas and form atmosphere s that will enha nce the cap - 
ture radius ( llnaba & Ikomal 120031: l^ nigawa & Oh tsukill20Tol) . 
According to the results of llnaba & Ikoma (i2003i) this will per- 



haps become important when oligarchs reach 0.1 MeQ Our ex- 
pression for the impact radius and collision rates, therefore, are 
lower limits when protoplanets are surrounded by a thick atmo- 
sphere. 

In summary, most of the mentioned effects become relevant 
for evolved (gas) planets only. Most damaging to our analysis 
are the pressure fluctuations that could virtually shut off accre- 
tion, or make it very difficult to model it analytically. However, 
this effect may only become effective for large planet masses 
(Eq. (l38ll). For lower mass planets, the build-up of a dense at- 
mospheres will enhance the accretion rates with respect to our 
prescription. We believe that for dynamically-cold protoplanets 
below 0.1 Me our prescription should be quantitatively correct if 
collective effects can be neglected. In future work we intend to 
test the validity of our analytic expressions in more convoluted 
environments that incorporate some of the above processes. 

7.2. Summary 

We have developed a framework for the calculations of particle- 
protoplanet interactions in a gaseous environment. This involves 
the integration of the equations of motions in the circularly re- 
stricted three body problem including drag forces. Using the 
above mentioned simplifications - most notably the assumption 
of a linear drag force, a smooth background density and head- 
wind velocity, and a 2D setting - we were able to reduce the 
general problem to a state that includes only two dimensionless 
parameters: the dimensionless headwind velocity (^i and the di- 
mensionless stopping time (Stokes number, St), see Sect. |2] A 
large parameter study of particle trajectories has been conducted, 
from which, as function of (uj and St, the impact radii are derived. 
We find that three accretion modes can be distinguished: 

- SettUng encounters. Particles settle to the protoplanet and the 
impact radius is independent of the size of the latter; 

- Hyperbolic encounters. Accretion proceeds like the usual 
gravitational focusing with the approach velocity being in- 
fluenced by gas drag. 

- (Drag enhanced) three body encounters. Interactions take 
place on the scale of the Hill radius and gas drag causes a 
fraction of the particles to become captured, which settle to 
the protoplanet. 

We have developed analytic recipes for all three encounters and 
found them to match very well to the results of our numerical 
study (except perhaps for the three-body regime). The recipe 
is summarized in Tableland Sect. 15.4! In Sect. |6] we have ex- 
tended our approach to the usual 3D-setting and calculated the 
accretion times of (proto)planets by sweepup of particles. We 
found that small particles are very unlikely candidates to grow to 
planetesimals of 10 - 10^ km, since their trajectories are tightly 
coupled to the gas. However, if the protoplanet has reached a 
size of ~10'' km it can very quickly accrete St ~ 1 particles 
through the settling mechanism. Since these fragments do not 
suffer from the protoplanet gravitational scattering (they quickly 
circularize), accretion under such conditions represents an av- 
enue for quick growth, especially in the outer disk - provided 
they are not lost by radial drift. 
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Appendix A: Asymptotic limits of Eq. (22) 

In this appendix we consider the asymptotic limits of Eq. ([22li 
and sh ow the correspondence to the findings of iKarv et aTl 
(Il993h . 

Three limits of Eq. ( l22l i can be identified 



P - 



4ff„^./St 
3al 



ifSt« 1; 

if St «: 1 and OpSt/^ui «: 1; 

if a,, Si/ (^, » 1. 



(A.l) 



These three regimes correspond to the cases where the square- 
root term of Eq. (l22l i evaluates to l/2St, 1, and 3apSt'/4-(u„ re- 
spectively. In all limits we have assumed that ap <sc (n,. 

Rewritten in physical units Eq. (lA.lb reads (i.e., we multiply 
by RhUh) 



pgeo 
col 



IRpVhsf, for St «: 1; 

ARpVh^ISi for St «; 1 and TJ^St/yh* 
for RpSilvh-^ » 1. 
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The interpretation of the first two limits is straightforward. If 
St «; 1 particles arrive with the headwind velocity, Whw In 
the 2D-setting the cross section is 2Rp, so the collision rate 
is IRpVhsf,- Similarly, in the second limit the particles approach 
from the x-direction at a speed of 2L)hw/St (see Eq. (I19ab ). 

In the third limit the particles approach once again from the 
j^-direction. However, for these very large particles, the approach 
velocity is given by the Keplerian shear (3Qx/2) instead of the 
headwind, see Eq. (I19bl i. The approach velocity at the point of 
intersection, i.e., at x - Rp, is then Va = 3RpQ./2. Multiplied 
by 2Rp this reduces to the given expression. The dependence on 
Rp may seem counter intuitive but is natural in situations that 
involve shear. 

The study of iKarv et al.l (1993) concerned massive particles 
(i.e., the third limit of Eq. ( IA.2l i). Kary et al. (1993) gave an an- 
alytic expression for the impact probability or efficiency Pgff of 
a particles while crossing the semi-major axis of the protoplanet 
due to radial drifto 



Peff = 



3/?2q2 
ISnKauh^ 



(A.3) 



where K is the drag constant for a drag law that quadratically 
depends on velocity. 

In our case Pgff can be found by taking the ratio of the colli- 
sion rate P^oi to the mass inflow rate \27Tav^\ of the particles. For 
the third limit of Eq. ( IA.2b we obtain 



pgeo 



col _ P 



InaVx A7iav\^y^ 



(A.4) 



where we used t hat Vx - 2i\yjlSl i n this regime. Equation (I A. 4b 
is different from lKarv et al.l (Il993h 's result for three reasons; 

1 . iKarv et all ( Il993h considers a drag force quadratic in veloc- 
ity where particles move at a drift velocity v^ - -2Ku^^/Q. 
(K has units cm"'). However, we can mimic the linear drift 
law by substitution of K = O/uhwSt (cf. Eq. (I19ab ) into 
Eq. (IX31 l.. 



'^ N ote that we give the dimensional form. Equation (7) of lKarv et al.l 
l ll993h is expressed in dimensionless units (but not in Hill units). 
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2. In our approach we have not accounted for the variation of 
the approach velocity over the impact range. More correctly, 
the mean approac h velocity is~i = 'iQ.RplA. 

3. iKarv et aTl (Il993h do not take account of the planetesimals 
coming from the negative i/-direction (the third quadrant), as 
they (correctly) argue that any such body should already have 
impacted during its approach from the first quadrant. This 
is due to the highly symmetrical setting in this limit, which 
our naive reasoning above does not account for. Equation 
Eq. ( IA.4l i is therefore too large by a factor of two. 



With these corrections Eqs. lA.3l and lA. 41 agree 



Appendix B: The settling path 

In the limit of St «; 1 we can use the approximation that the 
particle is always in the settling regime, i.e., its velocity is given 
by V = ^gts (in Hill units): 

y.^-^St, (B.la) 

Vy^-^Si-U (B.lb) 

r 

Thus, the particle path obeys the differential equation 

dy _vy _y ^ Co, jx^ + y^f'^ 
dx Vx X 3St X 

which is slightly simplified if expressed in angular coordinates 
by the substitution y = xtanO. Then, dy/dx = tanfl + x{l + 
tan^ 0)d0/dx and Eq. ( IB. 21 ) becomes in terms of x and 

d0 t 
x(l + tan- 0)— = ;^(1 + tan^ 0f^-, (B.3) 

dx 3St 



which is equivalent to 



d0 llu, 
COS0— = — X. (B.4) 

dx 3St 

Straightforward integration gives the solution 

sin6i= ;^;cHC, (B.5) 

oSt 

with C the integration constant, which we obtain by the require- 
ment that at = 7r/2 the particle starts out ■a.Vx — xq. Thus, the 
full solution to the orbit under these conditions is 

l-sin0=^(x2-x-). (B.6) 

During the encounter decrease from 7r/2 to -7r/2 and 1 - sin 
increases from to 2. The particle's x coordinate then decreases 
by an amount that depends on the numerical value of ^((,/6St and 
xq. The particle settles to the origin if x = can be reached. 
The largest value of xq for which this is possible is ^ = xo - 
-y/12St/^u,. This is the impact parameter ha-- 

For example, for St - 10"^ and ^„, - 1, we obtain ha - 
0.35, a value that is reasonably close to the numerically derived 
xj — 0.38 (see Fig.[5]l- For lower Stokes number the agreement 
becomes better. 

In Sect. 15.1.21 we have used the expression Au - 3St/h^ for 
the impulse change in the settling regime. For h - ^JT2Stf^ 
this corresponds to a velocity change of Av = fu,/4. We further 
argued that for settling encounter the approach velocity Ua should 
be changed by an amount Au ~ Va- For small Stokes values Va - 
(ui- Therefore, for St «; 1 the criterion for accretion by settling 
becomes Vg > ^wl^- In Sect. l5.1.2l we have applied this criterion 
generally. 



